function optimized_rates = find_unfolding_parameters_main_multistart(transitions_matrix, force_rate, force_distribution, detect_fluxes, varying_parameters,  initial_rates, lower_bounds, upper_bounds, chi2_minimization, min_method)

% input - transitions_matrix, NxN (N - total number of protein states
% during unfolding process, 0 in this matrix such transition does not
% exist, 1 - mark existing transitions. Index of the initial state = 1,
% protein is in this state prior to unfolding.\

%force_rate - in pN/s

%force_distribution - matrix, 2xM, data for unfolding force histogram;
%(1,:) - x axis [force], (2,:) - y axis [number of events]

%detect_fluxes - which fluxes contribute to unfolding force histogram

%varying_parameters - indexes of parameters which will be optimized

    if(size(transitions_matrix,1)==size(transitions_matrix,2) && length(varying_parameters)==length(lower_bounds) && length(varying_parameters)==length(upper_bounds))

        for i=1:size(transitions_matrix,1)
            transitions_matrix(i,i) = 0;
        end
        transitions_number = sum(sum(transitions_matrix));
        disp(['Totally ' num2str(transitions_number) ' transitions']);

        if(transitions_number==length(initial_rates)/2)
            f = @(x)minimization_function_v4(x, transitions_matrix, force_rate, force_distribution, detect_fluxes,  varying_parameters,  initial_rates, chi2_minimization);
   
           % i -motif:
%             initial_rates = [0.00010512   0.00002975  0.00017538      2.8725      1.6445      1.3095];     

            if(min_method==1)
                options = optimset('Display', 'iter', 'MaxIter', 1, 'TolFun', 0.01);
                [optimized_rates, fval] = fminsearch(f, initial_rates(varying_parameters), options)
            else
                problem = createOptimProblem('fmincon', 'objective', f, 'x0', initial_rates(varying_parameters), 'lb', lower_bounds, 'ub', upper_bounds);
                problem.options=optimoptions(problem.options,'UseParallel',true);
                ms = MultiStart('Display', 'iter', 'MaxTime', 56000,'UseParallel',true);
                [xg fg flg og] = run(ms, problem,1000)
                %gs = GlobalSearch('Display', 'iter', 'MaxTime', 56000);
                %[xg fg flg og] = run(gs, problem)
                optimized_rates = xg
            end
        
        else
            disp(['Transition matrix is wrong']);
        end
    else
        disp(['Transition matrix is not square or the length of lower_bounds/upper_bounds array is not equal to the length of varying_parameters array']);    
    end

end